# delimit ;

drop _all;
clear;
clear matrix;
clear mata;
set more off;
set matsize 11000;
set maxvar 11000;

* Citation: Rangel, M. A. and T. Vogl 
"Agricultural Fires and Health at Birth"
The Review of Economics and Statistics;

******************************************************************************
ABSTRACT

Do file employs birth outcomes, fetal deaths, death records, 
fires and weather measurement.

OFFICIAL SOURCES OF DATA:
INPE
CETESB
DATASUS

WE PROVIDE DATA IN STATA FORMAT:
fires_birthoutcomes_wfetaldeath.dta

******************************************************************************;

global data_repository "../RESTAT_data";
global logs_final "../RESTAT_logs";
global dofiles "../RESTAT_do";
local bootruns "249";

************************************************************************
Table 4 - causal inference fires and birth outcome + fetal deaths
************************************************************************;
 
use "$data_repository/fires_birthoutcomes_wfetaldeath.dta", clear;


*** schematics for coefficient sums;


foreach m in nr_risks50_week pv_nr_risks50_week res_nr_risks50_week
			 znr_risks50_week zpv_nr_risks50_week zres_nr_risks50_week			 {;

* first the coef on 0 is 13 weeks long; 
forvalues j=1(1)12 {;
replace `m'_`j'=`m'_`j'-`m'0;
};
* now do the same for all others so that you have a rolling window;
forvalues j=2(1)13 {;
replace `m'_`j'=`m'_`j'-`m'_1;
};
forvalues j=3(1)14 {;
replace `m'_`j'=`m'_`j'-`m'_2;
};
forvalues j=4(1)15 {;
replace `m'_`j'=`m'_`j'-`m'_3;
};
forvalues j=5(1)16 {;
replace `m'_`j'=`m'_`j'-`m'_4;
};
forvalues j=6(1)17 {;
replace `m'_`j'=`m'_`j'-`m'_5;
};
forvalues j=7(1)18 {;
replace `m'_`j'=`m'_`j'-`m'_6;
};
forvalues j=8(1)19 {;
replace `m'_`j'=`m'_`j'-`m'_7;
};
forvalues j=9(1)20 {;
replace `m'_`j'=`m'_`j'-`m'_8;
};
forvalues j=10(1)21 {;
replace `m'_`j'=`m'_`j'-`m'_9;
};
forvalues j=11(1)22 {;
replace `m'_`j'=`m'_`j'-`m'_10;
};
forvalues j=12(1)23 {;
replace `m'_`j'=`m'_`j'-`m'_11;
};
forvalues j=13(1)24 {;
replace `m'_`j'=`m'_`j'-`m'_12;
};
forvalues j=14(1)25 {;
replace `m'_`j'=`m'_`j'-`m'_13;
};
forvalues j=15(1)26 {;
replace `m'_`j'=`m'_`j'-`m'_14;
};
forvalues j=16(1)27 {;
replace `m'_`j'=`m'_`j'-`m'_15;
};
forvalues j=17(1)28 {;
replace `m'_`j'=`m'_`j'-`m'_16;
};
forvalues j=18(1)29 {;
replace `m'_`j'=`m'_`j'-`m'_17;
};
forvalues j=19(1)30 {;
replace `m'_`j'=`m'_`j'-`m'_18;
};
forvalues j=20(1)31 {;
replace `m'_`j'=`m'_`j'-`m'_19;
};
forvalues j=21(1)32 {;
replace `m'_`j'=`m'_`j'-`m'_20;
};
forvalues j=22(1)33 {;
replace `m'_`j'=`m'_`j'-`m'_21;
};
forvalues j=23(1)34 {;
replace `m'_`j'=`m'_`j'-`m'_22;
};
forvalues j=24(1)35 {;
replace `m'_`j'=`m'_`j'-`m'_23;
};
forvalues j=25(1)36 {;
replace `m'_`j'=`m'_`j'-`m'_24;
};
forvalues j=26(1)37 {;
replace `m'_`j'=`m'_`j'-`m'_25;
};
forvalues j=27(1)38 {;
replace `m'_`j'=`m'_`j'-`m'_26;
};
};



*** base model estimations;


sum meanfetaldeath lnlivepop [aweight=pop] ;
sum rate1day rate1week nr_all_bd1d   nr_fetalgrowth_bd1d  meanbr_apgar1 meanbr_apgar5   [aweight=poplive] ;

foreach depvar in  meanfetaldeath lnlivepop   {;

areg  `depvar'  
LOC*

NNE_wind_week0-NNE_wind_week_38
ENE_wind_week0-ENE_wind_week_38
ESE_wind_week0-ESE_wind_week_38
SSE_wind_week0-SSE_wind_week_38
SSW_wind_week0-SSW_wind_week_38
WSW_wind_week0-WSW_wind_week_38
WNW_wind_week0-WNW_wind_week_38
NNW_wind_week0-NNW_wind_week_38

nr_fires5_week0-nr_fires5_week_38
RH_dayavg_week0-RH_dayavg_week_38
TEMP_dayavg_week0-TEMP_dayavg_week_38
AVT_*
nonprvwinds_week0- nonprvwinds_week_38
nowindmeas_week0-nowindmeas_week_38

zpv_nr_risks50_week0-zpv_nr_risks50_week_38 
znr_risks50_week0-znr_risks50_week_38

[aweight=pop], absorb(date_BR) cluster(codestation);

outreg2 using "$logs_final/Table4_birthoutcomes13locales_`depvar'.out",  noaster replace;

};


foreach depvar in  rate1day rate1week nr_all_bd1d nr_fetalgrowth_bd1d meanbr_apgar1 meanbr_apgar5   {;

areg  `depvar'  
LOC*

NNE_wind_week0-NNE_wind_week_38
ENE_wind_week0-ENE_wind_week_38
ESE_wind_week0-ESE_wind_week_38
SSE_wind_week0-SSE_wind_week_38
SSW_wind_week0-SSW_wind_week_38
WSW_wind_week0-WSW_wind_week_38
WNW_wind_week0-WNW_wind_week_38
NNW_wind_week0-NNW_wind_week_38

nr_fires5_week0-nr_fires5_week_38
RH_dayavg_week0-RH_dayavg_week_38
TEMP_dayavg_week0-TEMP_dayavg_week_38
AVT_*
nonprvwinds_week0- nonprvwinds_week_38
nowindmeas_week0-nowindmeas_week_38

zpv_nr_risks50_week0-zpv_nr_risks50_week_38 
znr_risks50_week0-znr_risks50_week_38

[aweight=poplive], absorb(date_BR) cluster(codestation);

outreg2 using "$logs_final/Table4_birthoutcomes13locales_`depvar'.out",  noaster replace;

};


	*** Frisch-Lovell wild-bootstrap;	
	
	foreach depvar in  zpv_nr_risks50_week0 zpv_nr_risks50_week_13 zpv_nr_risks50_week_26  meanfetaldeath lnlivepop  {;

		quietly areg  `depvar'  
			LOC*

			NNE_wind_week0-NNE_wind_week_38
			ENE_wind_week0-ENE_wind_week_38
			ESE_wind_week0-ESE_wind_week_38
			SSE_wind_week0-SSE_wind_week_38
			SSW_wind_week0-SSW_wind_week_38
			WSW_wind_week0-WSW_wind_week_38
			WNW_wind_week0-WNW_wind_week_38
			NNW_wind_week0-NNW_wind_week_38

			nr_fires5_week0-nr_fires5_week_38
			RH_dayavg_week0-RH_dayavg_week_38
			TEMP_dayavg_week0-TEMP_dayavg_week_38
			AVT_*
			nonprvwinds_week0- nonprvwinds_week_38
			nowindmeas_week0-nowindmeas_week_38

			
			zpv_nr_risks50_week_1-zpv_nr_risks50_week_12 zpv_nr_risks50_week_14-zpv_nr_risks50_week_25 zpv_nr_risks50_week_27-zpv_nr_risks50_week_38
			znr_risks50_week_1-znr_risks50_week_12 znr_risks50_week_14-znr_risks50_week_25 znr_risks50_week_27-znr_risks50_week_38

			znr_risks50_week0 znr_risks50_week_13 znr_risks50_week_26
			
			[aweight=pop], absorb(date_BR) cluster(codestation);
			predict res_`depvar', residuals;
	
	};	
	
		foreach depvar in   meanfetaldeath lnlivepop {;

		preserve;
		
		log using "$logs_final/Table4_birthoutcomes13locales_`depvar'.log", replace;
	
	    set seed 20160817;
	
		adopath + "$dofiles";
			
	
		wild_areg res_`depvar'  
		res_zpv_nr_risks50_week0 res_zpv_nr_risks50_week_13 res_zpv_nr_risks50_week_26
		[aweight=pop]
		,ivar(codestation) tvar(date_BR) 
		bootstraps(`bootruns');

		
		
		adopath -  "$dofiles";
		
		drop res_`depvar';
		
		log close;
		
		restore;

	};


	drop res_zpv_nr_risks50_week0 res_zpv_nr_risks50_week_13 res_zpv_nr_risks50_week_26;
	
	
	
	foreach depvar in  zpv_nr_risks50_week0 zpv_nr_risks50_week_13 zpv_nr_risks50_week_26  rate1day rate1week nr_all_bd1d nr_fetalgrowth_bd1d meanbr_apgar1 meanbr_apgar5  {;

		quietly areg  `depvar'  
			LOC*

			NNE_wind_week0-NNE_wind_week_38
			ENE_wind_week0-ENE_wind_week_38
			ESE_wind_week0-ESE_wind_week_38
			SSE_wind_week0-SSE_wind_week_38
			SSW_wind_week0-SSW_wind_week_38
			WSW_wind_week0-WSW_wind_week_38
			WNW_wind_week0-WNW_wind_week_38
			NNW_wind_week0-NNW_wind_week_38

			nr_fires5_week0-nr_fires5_week_38
			RH_dayavg_week0-RH_dayavg_week_38
			TEMP_dayavg_week0-TEMP_dayavg_week_38
			AVT_*
			nonprvwinds_week0- nonprvwinds_week_38
			nowindmeas_week0-nowindmeas_week_38

			
			zpv_nr_risks50_week_1-zpv_nr_risks50_week_12 zpv_nr_risks50_week_14-zpv_nr_risks50_week_25 zpv_nr_risks50_week_27-zpv_nr_risks50_week_38
			znr_risks50_week_1-znr_risks50_week_12 znr_risks50_week_14-znr_risks50_week_25 znr_risks50_week_27-znr_risks50_week_38

			znr_risks50_week0 znr_risks50_week_13 znr_risks50_week_26
			
			[aweight=poplive], absorb(date_BR) cluster(codestation);
			predict res_`depvar', residuals;
	
	};	
	
		foreach depvar in    rate1day rate1week nr_all_bd1d nr_fetalgrowth_bd1d meanbr_apgar1 meanbr_apgar5 {;

		preserve;
		
		log using "$logs_final/Table4_birthoutcomes13locales_`depvar'.log", replace;
	
	    set seed 20160817;
	
		adopath + "$dofiles";
			
	
		wild_areg res_`depvar'  
		res_zpv_nr_risks50_week0 res_zpv_nr_risks50_week_13 res_zpv_nr_risks50_week_26
		[aweight=poplive]
		,ivar(codestation) tvar(date_BR) 
		bootstraps(`bootruns');

		
		
		adopath -  "$dofiles";
		
		drop res_`depvar';
		
		log close;
		
		restore;

	};	
exit; 

































*/

*** wild bootstrap routines;


foreach depvar in  meanfetaldeath lnlivepop {;

	preserve;
		
		log using "$logs_final/Table4_birthoutcomes13locales_`depvar'.log", replace;
	
	    set seed 20160817;
	
		adopath + "$dofiles";
				
		* display "`depvar'";
	
		wild_areg `depvar'  
		zpv_nr_risks50_week0 zpv_nr_risks50_week_13 zpv_nr_risks50_week_26 znr_risks50_week0 znr_risks50_week_13 znr_risks50_week_26 
		[aweight=pop]
		,ivar(codestation) tvar(date_BR) 
		control(LOC* zpv_nr_risks50_week_1-zpv_nr_risks50_week_12 zpv_nr_risks50_week_14-zpv_nr_risks50_week_25 zpv_nr_risks50_week_27-zpv_nr_risks50_week_38
		nr_risks50_week_1-nr_risks50_week_12 nr_risks50_week_14-nr_risks50_week_25 nr_risks50_week_27-nr_risks50_week_38
		NNE_wind_week0-NNE_wind_week_38
		ENE_wind_week0-ENE_wind_week_38
		ESE_wind_week0-ESE_wind_week_38
		SSE_wind_week0-SSE_wind_week_38
		SSW_wind_week0-SSW_wind_week_38
		WSW_wind_week0-WSW_wind_week_38
		WNW_wind_week0-WNW_wind_week_38
		NNW_wind_week0-NNW_wind_week_38

		nr_fires5_week0-nr_fires5_week_38
		RH_dayavg_week0-RH_dayavg_week_38
		TEMP_dayavg_week0-TEMP_dayavg_week_38
		AVT_*
		nonprvwinds_week0- nonprvwinds_week_38
		nowindmeas_week0-nowindmeas_week_38) 
		bootstraps(`bootruns');
				
		log close;
		
		adopath -  "$dofiles";
		
		restore;

};		



foreach depvar in rate1day rate1week nr_all_bd1d nr_fetalgrowth_bd1d meanbr_apgar1 meanbr_apgar5 {;

	preserve;
		
		log using "$logs_final/Table4_birthoutcomes13locales_`depvar'.log", replace;
	
	    set seed 20160817;
	
		adopath + "$dofiles";
				
		* display "`depvar'";
	
		wild_areg `depvar'  
		zpv_nr_risks50_week0 zpv_nr_risks50_week_13 zpv_nr_risks50_week_26 znr_risks50_week0 znr_risks50_week_13 znr_risks50_week_26 
		[aweight=poplive]
		,ivar(codestation) tvar(date_BR) 
		control(LOC* zpv_nr_risks50_week_1-zpv_nr_risks50_week_12 zpv_nr_risks50_week_14-zpv_nr_risks50_week_25 zpv_nr_risks50_week_27-zpv_nr_risks50_week_38
		nr_risks50_week_1-nr_risks50_week_12 nr_risks50_week_14-nr_risks50_week_25 nr_risks50_week_27-nr_risks50_week_38
		NNE_wind_week0-NNE_wind_week_38
		ENE_wind_week0-ENE_wind_week_38
		ESE_wind_week0-ESE_wind_week_38
		SSE_wind_week0-SSE_wind_week_38
		SSW_wind_week0-SSW_wind_week_38
		WSW_wind_week0-WSW_wind_week_38
		WNW_wind_week0-WNW_wind_week_38
		NNW_wind_week0-NNW_wind_week_38

		nr_fires5_week0-nr_fires5_week_38
		RH_dayavg_week0-RH_dayavg_week_38
		TEMP_dayavg_week0-TEMP_dayavg_week_38
		AVT_*
		nonprvwinds_week0- nonprvwinds_week_38
		nowindmeas_week0-nowindmeas_week_38) 
		bootstraps(`bootruns');
				
		log close;
		
		adopath -  "$dofiles";
		
		restore;

};		


